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I. INTRODUCTION AND MOTIVATION 

In hadronic collisions at ultra-relativistic energies the final state is quite involved in terms 
of the type, the number and the distribution of the produced particles. The extraction of the 
dominant physical mechanisms in such processes is a demanding task and in order to achieve the 
best possible understanding it is necessary to study many observables in wide kinematic regimes. 
For instance, significant attention has been given to collisions between light and heavy hadrons, 
like deuteron-gold at RHIC and the forthcoming proton-lead at the LHC and in both cases two of 
the most representative observables are related to single and double inclusive particle production. 
Considering the single inclusive particle production in the deuteron fragmentation region it has 
been observed, already a few years ago at RHIC [HE], a suppression compared to the cross section 
obtained by taking the nucleus as an incoherent superposition of A 1 / 3 nucleons. More recently, and 
regarding the double inclusive particle production, it has been seen an increasing suppression of 
the azimuthal correlation of the two hadrons when their transverse momenta are a few GeV [3J S] , 
as we move towards the deuteron fragmentation region. 

This is precisely the kinematic regime which encourages the search for parton saturation in 
the wave-function of the heavy hadron, the large nucleus. Natural qualitative and quantitative 
descriptions of the RHIC data and predictions for the LHC upcoming ones based on such a physical 
mechanism, along with the corresponding formulations, already exist for the nuclear modification 
factor R p a [SHE], which is related to the single inclusive cross section. Similarly, the di-hadron 
azimuthal correlations at RHIC offer a unique environment to test parton saturation [TT} [T8] and 
in fact the corresponding data have been understood in that context [T9H22] . 

Thus, one is particularly interested in hA — >► h'X and hA — >► h\h2X, with h a projectile 
hadron whose wavefunction is not saturated, like a proton at not too high energy so that its 
small- x evolution can be neglected, and A a target who can be dense, like an ultra-relativistic 
heavy nucleus. Let us again look at single particle production first for which the corresponding 
diagram at the partonic level, say for quark production, is shown in Fig. [lj(a). A large- x quark 
from the projectile interacts via multiple gluon exchanges with the small- x components of the 
target and then it is measured in the region which is forward in (pseudo) rapidity. The target is 
viewed as a potentially large color field the Color Glass Condensate (CGC) (see e.g. [23J) and 
the interaction of a parton with transverse position x with such a field is described by a Wilson 
line along its trajectory. Taking the modulus squared of the diagram [I] (a) in coordinate space, 
averaging over initial colors and summing over final ones, we find the cross section qA — » qX to 
be given by the Fourier transform of a color dipole 5, that is, a trace of two Wilson lines in the 
fundamental representation, which is an overall colorless object. 

The above discussion naturally extends to the case of double particle production when both 
particles are detected at the same rapidities. Without any loss of generality, les us focus on qg 
production^] at the partonic level with the respective diagrams shown in Fig. [11(b) and (c). The 
large- x projectile quark splits into a quark-gluon pair, and this splitting can occur either before 

1 Here and in the previous paragraph, we discuss only representative cases of single and double inclusive particle 
production which are taken to be quark and quark-gluon production respectively. Other possibilities like gluon 
and gluon-pair production have been also studied and understood 17, 24-28 . 
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or after the interaction with the target. Now one arrives at an inclusive qA — >> qgX cross section 
|18| which is given by a Fourier transform of various terms involving two, four or six Wilson lines. 
Such a term with the maximal number of Wilson lines is SQ, that is, a color dipole times a color 
quadrupole, where the latter is a trace of four Wilson lines, and in fact it is not hard to understand 
the counting of Wilson lines. For instance, diagram (b) involves one in the fundamental and one in 
the adjoint representation and the latter can be expressed in terms of two fundamental ones. When 
multiplied with its complex conjugate it gives rise to the aforementioned SQ term which involves 
six Wilson lines (plus a 1/N% correction). Clearly, one realizes that the production of more and 
more particles at forward rapidities will involve more and more Wilson lines. Interestingly enough, 
it was recently shown that in the large- N c limit any forward multi-particle production cross section 
can be expressed only in terms of dipoles and quadrupoles [29] . 

Thus, in general, one needs to calculate correlators of the form (O)y, where O is constructed 
from such multipoles. Obviously the Wilson lines depend on the target field and the average has to 
be taken with the target probability distribution Wy[»A], which simply gives the probability to find 
a given configuration in the target wavefunction. The rapidity Y is determined by the kinematics 
of the process under consideration; for example, for qg production at forward rapidities one has 
x = exp(— Y) — [\k\ exp(— \yk\) + \q\ exp(— \y q \)]/yfs, with q and y q the transverse momentum and 
(pseudo) rapidity of the produced quark, k and those of the gluon and the center of mass 
energy. At moderately small values of x one typically invokes the McLerran-Venugopalan (MV) 
model [HOI El], which is equivalent to a Gaussian wavefunction Wy[a] and thus allows for explicit 
calculation of high-point correlators [TTJ [321 [33]. 

As we move towards higher values of y q and y^, we probe components of the target with smaller- 
x and it becomes necessary to take into account the evolution in Y of Wy[^4]. Even though this 
is a classical probability distribution its evolution is quantum and satisfies the JIMWLK equation 
[31H38] which arises from the resummation of the dominant terms in ln(l/x) in the presence of a 
potentially strong color field. Using the Langevin form of the JIMWLK equation [39], a critical 
work on its numerical solution was recently performed, where the very good accuracy of a particular 
Gaussian approximation, the extension to arbitrary Y of the MV model, was observed for some 
configurations of high-point correlators [40J. Immediately after it was analytically understood why 
the approximation scheme based on a Gaussian wavefunction provides a quasi-exact solution to 
the JIMWLK equation gUHZJ. 

Our effort, which is mostly numerical, may be considered as complementary to the work in [40] 
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to some extent, but it is completely independent for a number of reasons, (i) The method here is 
different; instead of studying JIMWLK as a Langevin equation on a lattice, we directly analyze 
evolution equations for given configurations. But even if one starts from a simple comfiguration, 
evolution always creates more general ones (cf. Figs. [8] and [9]) and eventually our tests probe a 
much wider sample, (ii) We shall explore a larger kinematical space and as a byproduct we shall 
be able to give a numerical verification of the Levin- Tuchin formula [43] (see also [441447] ). (iii) We 
will study both fixed and running coupling evolution, where the scheme in the latter case can be 
chosen at our will, (iv) We can consider any value for the number of colors 7V C , even though we 
shall mainly work in the multicolor limit. 

The paper is mostly devoted to technical aspects concerning the validity of a certain approx- 
imation scheme and its structure is the following: In Sect. [II] we review the JIMWLK equation 
and how it determines the evolution of multi-parton correlations. In this Section we also introduce 
most of our definitions along with the corresponding notation. In Sect. Ill we give an accurate 
numerical solution to the Balitsky-Kovchegov (BK) equation, that is, the large- N c equation for 
the dipole [IHJ S9]. In particular we focus on the regime where saturation has been reached, a 
necessary step for our purposes, and we verify for the first time the validity of the Levin- Tuchin 
formula. Moreover, from the form of the solution at saturation, we also confirm the two most 
dominant terms in the asymptotic expansion of the saturation momentum Q s for fixed coupling 
evolution. In Sect. [IV] we briefly reflect on the Gaussian approximation and why it is expected 
to work accurately [HJ H2]. In Sect. [V] we study numerically a Mean Field Equation arising from 
the Gaussian approximation to the JIMWLK Hamiltonian and show that it is in agreement with 



the BK equation to logarithmic accuracy at saturation as expected. In Sect. VI we examine two 



quadrupole configurations and show, for both fixed and running coupling evolution (the latter in 
the Balitsky, smallest dipole and daughter dipole prescriptions), that the accuracy is not restricted 
to be logarithmic, giving further evidence that the Gaussian approximation is a "quasi-exact" solu- 
tion. We point out a potential deficiency of a "too simple" Gaussian approximation (in particular 
the extrapolation to arbitrary Y of the MV model) in the running coupling case and at saturation, 
but it seems not to be crucial for practical purposes. In Sect. VII we outline how one can calculate 
correlators with open color indices and finally in Sect. VIII we conclude. 



II. THE JIMWLK EQUATION AND MULTI-PARTON CORRELATIONS 

The Color Glass Condensate (CGC) is a modern effective theory for the small- x components 
of the wavefunction of an ultra-relativistic hadron. It relies on the idea that gluons which carry a 
small fraction x of the hadron's longitudinal momentum can be described as a random distribution 
of classical color fields generated by sources with larger momentum fractions. As a result of the 
high energy kinematics, the distribution of the color sources is frozen due to Lorentz time dilation 
and the color field, in a suitable gauge, has a single non-zero component. More precisely, and using 
the standard definitions for the light-cone coordinates x^ = (x + ,x _ ,x), with x ± = (t±x 3 )/V2 
and x = (x 1 , x 2 ), this gauge field is A%(x) = 5 fi+ a a (x~ , x) if the hadron moves along the positive 
x 3 -axis. The CGC weight function Wy[a] gives the probability that the hadron be described by 
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the configuration a and it is a functional probability distribution whose knowledge allows the 
determination of the correlations of the gauge field. The latter contain all the detailed information 
about the evolution of the hadron with increasing rapidity Y = ln(l/x), from an initial value 
Yq = ln(l/xo) to the value Y of interest. At high energy where a s (Y — Yq) > 1 and to leading 
order with respect to the large logarithm ln(xo/x), this evolution obeys a renormalization group 
equation, the JIMWLK equation [341438] . In a Hamiltonian form it reads 

^W Y [a]=HW Y [a], (2.1) 

where H is the JIMWLK Hamiltonian and it is a second-order functional differential operator 
whose most elegant and convenient for our purposes form was given is [50] 

ab S 5 



Here we have used the economical notation f u = J d 2 u . . . , defined the dipole kernel M. [51] 



(u — v) 2 
{u — z) 2 (z — v) 



and introduced the Wilson lines 



vH = p 



exp 



\g J dx a a (x ,x)T a 



(2.4) 



where T a is in the adjoint representation and with P denoting path-ordering in x . The precise 



action of the functional derivatives appearing in Eq. (2.2) will be explained later on. The above 
form of the Hamiltonian is valid only when acting on gauge-invariant functionals of a a: like gauge- 
invariant products of Wilson lines. This will be the case for most of our analysis with exceptions 
to be discussed at the end in Sect. |VIl| Needless to say, in order to specify our problem completely, 
we need an initial condition for Eq. ( |2.1[ ) at Y = Yq. At least for a sufficiently large nucleus 
(A ^> 1, with A the atomic mass number), this initial condition is typically provided by the 
McLerran-Venugopalan (MV) model [30|. l3Tj . 

Physical observables are represented by gauge invariant operators 0[a] constructed with the 
gauge color field a a . Their expectation value can be computed as a functional average with the 
CGC weight function, that is 

(0) Y = J VaO[a]W Y [a}. (2.5) 

The above makes clear that, even though Wy [a] is obtained by a quantum calculation, the averaging 



procedure is classical. Differentiating the above with respect to y, using Eq. (2.1), and finally 
integrating twice by parts, we arrive at the evolution equation 

^ - (HO) Y (2.6) 

for the observable under consideration. This is not a functional equation anymore, but an integro- 
differential equation as we shall see in a while in specific examples. Still, this is not much easier 
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to deal with since, due to the non-linear dependence of the Hamiltonian (2.2) on the field a a , 



Eq. (2.6) is in general not a closed equation, but just a member of an infinite hierarchy of coupled 



equations, the Balitsky equations [IE1|52]. The JILWLK equation (2.1) and the Balitsky hierarchy 
offer complementary views on the high-energy evolution. On the one hand, Eq. (2.1) describes the 



evolution of the target by the emission of an additional gluon with rapidity between Y and Y + dy 
in the background of the gauge color field a built by previous emissions, at rapidities smaller than 



Y. The Wilson lines in the Hamiltonian (2.2) correspond to the propagation of this new gluon in 
the background field and in the eikonal approximation. On the other hand, the Balitsky hierarchy 
focuses on the projectile evolution and in particular on the operator describing its scattering off the 
target. This scattering is again considered in the eikonal approximation and therefore the operator 
O is naturally constructed with Wilson lines, where each one of them represents a parton in the 
projectile. 

We shall mostly focus on the color dipole made by a quark-antiquark pair in an overall color 
singlet state, with the ^-matrix 



S XlX2 — Sx}x 2 — t r (^aci V^)' 



(2.7) 



and the color quadrupole, a system of two quarks and two antiquarks also in a color singlet, for 
which 



Qx!X 2 X 3 X4 — SxiX 2 X 3 X4 ~ jy ^ T (Yxi Kc 2 ^X 3 K» 



X4 j 



(2-f 



Here W and V are Wilson lines like in Eq. (2.4), but in the fundamental representation. In general 



one can consider projectiles made with n quarks and n antiquarks, for which 
S XlX \...x 2n -iX2n — -jy - tr(V^ 1 T4. 2 . . . V X2n _ 1 V X2n ), 



(2.9) 



and eventually high energy evolution mixes such single-trace operators with multi-trace ones of the 
form 



£>=Y tr« ^ • • .)^tr(VlV y2 . . . )i- tv(VlV Z2 ...). 



N, 



(2.10) 



The quadrupole denned above in Eq. (2.8) is the first non-trivial operator which can probe multi- 



parton correlation in the target wavefunction and will serve as the prototype in our study. In order 



to construct evolution equations according to Eqs. (2.2) and (2.6), it is necessary to specify how 



the functional derivatives w.r.t. a a act on observables. They act on endpoints of the Wilson lines 
according to 



6a" 



V* = ig6 xu t a Vl 



6a" 



V x = -ig6 xu V x t a , 



(2.11) 



with t a in the fundamental representation and where we introduced the shorthand notation 5 XU = 
x — u). Because of their action on W (by convention), they are called "left" derivatives and 



could be also denoted as S/Sa^. Using these rules within Eqs. (2.2) and (2.6), it is straightforward 
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(a) (b) (c) (d) 

FIG. 2. (a) A real emission of a gluon from a color dipole. (b) A virtual emission, (c) and (d) The 
corresponding diagrams in the large- N c limit. In all diagrams the dashed line stands for the interaction with 
the target. 



to derive the evolution equations satisfied by the S'-matrices for the dipole and the quadrupole. 
The resulting equation for the dipole is 

9(S XiX2 )y 



dY 



a r 

L 



M 



X1X2Z 



ZX2 



Sx\X2)Y 1 



(2.12) 



where we have defined a = a s N c /7r. Even though derived by evolving the target, this equation has 
an easy interpretation in terms of projectile evolution, as shown in Fig. [2] The quadratic term in 
S has been generated by the real part of the Hamiltonian H reai \, that is, by the last two terms in 
the parenthesis in Eq. ( |2.2| ). It describes the splitting of the original dipole (x\ 1 X2) into two new 
dipoles (a?i, z) and (z, X2), which subsequently scatter off the target. More precisely, the evolution 
step consists of the emission of a soft gluon, hence the original dipole gets replaced by a quark- 
antiquark-gluon system, but in the large- N c limit this emission is equivalent to the aforementioned 
dipole splitting. The negative, linear in S, term has been produced by the virtual part of the 
Hamiltonian i? v irt 5 that is, by the first two terms in the parenthesis in Eq. (2.2) and corresponds 
to the reduction in the probability for the dipole to survive in its original state. Notice that color 
transparency requires S xx = 1 and thus the potential short-distance singularities, arising from the 
dipole kernel M, at x\ = z and X2 = z cancel between the real and virtual terms. 

A word of caution should follow here. Eq. (2.12) is valid for any value of N C1 but still, both 
terms on the right hand side are of the same order in N c . In fact, terms suppressed by 1/N% have 
been generated in intermediate steps of the calculation but they have canceled in the final result. 
More precisely, when acting with i7 rea i on S XlX2 we get a contribution 
1 a 
2^ 

which cancels with an opposite in sign contribution coming from the action of H v [ rt . 
Similarly one can derive the evolution equation for the quadrupole, which reads 

9 {Qx 1 X2 Xs X4 ) Y 



<M-X\X2zS XlX2 1 



(2.13) 



dY 



a 
4n 



(M X1X2Z + M 



■ r x,\z -M. X2X4Z )(S XlZ Q ZX2X3X4 

)y 



H~ (>M XlX2Z -\- J\ / i X2XsZ 

H" (>M-X2X3Z H~ -M-X^X/^Z 

+ A4. X3X4Z 



>MxiX3z){S Z X2QxiZX3X4 : )y 
*M-X2X4z)(Sx3zQx\X2ZX4)y 

-M- x 1 X3 z ) ( S ZX4 Q x 1 X2 X3 z 

)y 



(>M. XlX2Z + -M-X3X4Z ~\~ >M-xiX4Z ~\~ -M-x 2 X3z)(QxiX2X3X4 )y 
- (M Xl X 2 Z + M X3 X4Z ~ M 

~ (>M.xiX4Z + M XQ x 3 z 



IXIX3Z 
- M Xl X 3 Z - M 



M. X2X4Z )(S XlX2 S X3X4 

)y 



X2X4Z 



)(S X3X2 S xiX4 )y • (2.14) 
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(a) (b) 

FIG. 3. (a) A real emission of a gluon from a color quadrupole in the large- N c limit corresponding to the 
term proportional to M. X2X3Z S ZX2 Q XlZX3XA . (b) A virtual emission corresponding to the term proportional 
to —M Xl x 3Z S X3X2 S XlX4 . In both diagrams the dashed line denotes the interaction with the target. 



Even though this looks considerably more involved than the dipole equation, a similar discussion 
applies and two representative diagrams are shown in Fig.jiJ The terms involving (SQ)y m the right 
hand side are real terms describing the splitting of the original quadrupole into a new quadrupole 
plus a dipole, and have generated by the action H rea \. The virtual terms involving (Q)y an d (SS)y 
are necessary for probability conservation, and have been generated by H v [ Tt . Once again, all terms 
subleading at large N c , separately generated by the two parts of the Hamiltonian, have canceled 
in the final equation, and all the short-distance singularities of the dipole kernels at X{ — z, with 
i = 1, 2, 3, 4, cancel among the various terms. 

All the above features generalize to the evolution equations obeyed by the single-trace ob- 
servables given in Eq. (2.9). As already evident in Eqs. (2.12) and (2.14), these equations are 
generally not closed and they couple single-trace observables with multi-trace ones. For instance, 
the quadrupole equation involves the 4-point function (SS)y and the 6-point function (SQ)y, 
which in turn will couple to even higher-point correlators. The equations obeyed by the multi- 
trace observables exhibit the new feature that they involve explicit 1/N% corrections. These arise 
when two functional derivatives in Eq. ( |2.2[ ) act on Wilson lines which belong to different traces 
(see e.g. Appendix F in [53J for an example). At large N c , one can neglect these terms and check 
that the hierarchy admits the factorized solution 

(d)y * <± trOtf, V„ . . . «(VlV K . . . )) r (±*iKV- ■••))„..•■ (2-15) 

so long as this factorization is already present in the initial condition. Therefore the hierarchy 
simplifies in a drastic way as it decomposes into a set of equations which can be solved, in principle, 
one after the other. More precisely, Eq. (2.12) becomes a closed equation for (S)y, the well-known 
BK equation [IE! 09], Eq. (2.14) becomes an inhomogeneous equation for (Q)y whose coefficients 
depend on (S)y [17J, and so on. Good analytical understanding and reliable numerical solutions 
to the BK equation exist, and we will discuss these matters in the next section. However, already 
starting from the quadrupole, it seems difficult to numerically deal with the large number of 
variables in the higher-point equations and the non-locality in the transverse space at the same 
time. In that case one can rely on an alternative formulation of the JIMWLK evolution as a 
Langevin equation [39] which can be solved on a lattice [IUJ |5U [55], or develop well- motivated 
approximate schemes and we shall turn our attention to the latter in Sect. |IV| 
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III. SOLUTION TO THE BK EQUATION AND THE LEVIN-TUCHIN FORMULA 



Let us briefly review the analytic solution to the BK equation. We shall assume a homogeneous 
target so that the amplitude (T XlX2 )y = 1 — (S XiX2 )y depends only on the magnitude r = T\2 — 
\x\ — #2 1- Still, the solution is not analytically known, but one can construct a piecewise one 
when aY > 1, by considering two regimes. The regime r < 1/Q S where the target is dilute and 
the scattering weak and the regime r > 1/Q S where the target is dense and the scattering strong. 
The borderline in between corresponds to the saturation momentum Q s and can be determined 
from the solution to the BFKL equation, that is, the linearized in T version of the BK equation, 
supplemented by appropriate boundary conditions. For fixed coupling one finds that the energy 
dependence of the saturation momentum is determined by [56l [57] 
1 dlnQ 2 = X ( 7g ) 3 i 
2 7s aY' 



a 



dY 7s 

where the "anomalous dimension" 7 S related to saturation is determined by 

X'ils) = Xils)hs => Is ~ 0.628. 
In the above, x{l) IS the eigenvalue function of the BFKL equation [53 E0] given by 

X( 7 ) = 2^(1) - V(7) - ^(1 - 7), 



(3.1) 



(3.2) 



(3.3) 



with ip the logarithmic derivative of the T-function. The amplitude on this side of the saturation 
line, that is for r <C 1/Q 5 , reads [56| [57] 



(f) y = (r 2 Q2)7 Mn 



r 2 Q 



+ c exp 



lnVQ^ 
DmY 



(3.4) 



an expression which is valid in the region < 1/r 2 < Q 2 S ex.p(D s aY), where c is a constant of 
order (9(1) and D s = 2x // (t 5 ) « 97 is the diffusion coefficient. When | ln(r 2 Q^)| < \/D s oY the last 



factor in Eq. (3.4), which describes diffusion, can be set equal to unity and the amplitude exhibits 
geometrical scaling [56 l [57 f I6TH63] ; it depends only on the combined variable r 2 Q 2 s . 

Let us now look at what happens when r ^> 1/Q S . In this regime the 5-matrix approaches 
its black-disk limit, i.e. (S)y 0, and thus we can neglect the term quadratic in (S)y m the 
BK equation. To do this properly, we need to restrict the region of integration in the transverse 
coordinates to l/Q 2 s <C \x{ — z\ 2 <C r 2 , with i = 1,2. The lower limit emerges as the boundary 
determining the transition from a region where the scattering is weak to a region where it becomes 
strong. The upper limit determines the dominant logarithmic contribution to the r.h.s. of the BK 
equation. We easily find that 



d(S) 



Y 



dY 



—a 



I/O: 



dz 2 - 



-aln(r 2 Q 2 s )(S) Y , 



2 Z 



(3.5) 



and using the leading behavior for the rapidity dependence of Q s given in Eq. (3.1) we convert the 
derivative w.r.t. Y to one w.r.t. \u(r 2 Ql). Then it becomes trivial to solve Eq. (3.5) and we arrive 
at 



(S)y ~ exp 



Is 



Mis 



ln 2 (r 2 Q 



(3.6) 
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FIG. 4. Numerical extraction of the coefficients in the asymptotic expansion (3.9) at saturation and for 
fixed coupling. The coefficients are plotted as a function of the dipole size r. 

where it should be reminded that only this leading term in the exponent is under good control. 



Eq. (3.6) is equivalent in its functional form to the Levin- Tuchin formula [43] (see also 
apart from the coefficient in the exponent which is different. Let us also add here, that this 
coefficient is further modified when N c is finite, namely the exponent should be multiplied by a 



factor 2Cf/N c = (N 2 — 1)/^? ESSE]. Again, geometric scaling is manifest in Eq. (3.6) and in 
fact this property is valid everywhere in the region Aq CD <l/r 2 <Q^ exp[y/D s aY]. Even though 
we do not have an analytical solution to cover the whole region, one can construct appropriate 



interpolations. For instance, and neglecting also the second factor in Eq. (3.4) which has only a 
weak logarithmic dependence, such a convenient expression is 

(S) Y = exp {-— ^-y In 2 [1 + b(rQ a y] } , (3.7) 

with b a constant of order 0(1). 



Going back to Eq. ( |3.6| ) and using Eq. (3.1), let us note that, for fixed r, the dominant term in 
the exponent is proportional to the square of Y. This has a direct physical interpretation which 
should have been obvious from the derivation. One factor of Y already appears in zero-dimensional 
particle models where all transverse coordinates are suppressed and represents the fact that dP/dY 
is proportional to — P, where P is the probability that a particle does not split. An extra factor 
of Y arises in QCD because the available phase space for an emission of a gluon increases with 
HQl) ~ Y. 



Here we shall put forward the task to numerically verify Eq. (3.6) for two reasons. The first is 
simply that such a study has not been done so far (see [44J for a related work). The second is that 
it is necessary for the comparisons to be performed in Sect.[VJ To this end, by also using Eq. (3.1), 
we can write the large- y and large-ln(l/r 2 ) expansion 

We compare this expression with fits of the rapidity dependence of numerical solutions of the BK 
equation^] with a function of the form 

- \n(S) Y = ci (aY) 2 + c 2 (aY) + c 3 + c 4 (aY) ln(aY). (3.9) 

2 Details about the numerical implementation of the BK equation and the other evolution equations considered later 
in this article are given in appendix. 
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Coefficient 




c 2 /ln(l/r 2 ) 


c 3 /ln 2 (l/r 2 ) 


c 4 


Analytical 


X ( 7s )/27 S = 2.442 


-1 


7,/2x(7 s ) = 0.1024 


-3/2 7s = -2.390 


Numerical (5 < Y < 20) 


2.466 


-1.002 


0.1029 


-2.601 


Numerical (20 < Y < 200) 


2.443 


-1.000 


0.1045 


-2.313 



TABLE I. Analytical and numerical values of the coefficients in the asymptotic expansion of — ln(5(r))y 
at saturation as in Fig. [4j The agreement of the numerics in Fig. [4] with Eq. (3.8) is excellent. 



The excellent numerical verification of all terms in (3.8) is demonstrated in Fig. [4] and in Table [IJ 



This makes clear that Eq. (3.6), including the exact coefficient in the exponent, gives the correct 
approach to the unitarity limit within the context of the BK equation. Notice that the asymptotics 
sets quite fast and the Levin- Tuchin law is valid also for "reasonable" values of Y, albeit this 
happens for rather large dipoles. Furthermore, it is remarkable that, without focusing in the 
transition region around Q s , we get a numerical confirmation for the first two terms of its energy 



dependence as given in Eq. (3.1). 

It is not hard to see how the above description is modified when we consider the running coupling 
BK equation [64-66J. In general, the coupling runs according to 
1 . . i t UN c -2N f 



a(r) 



with b 



(3.10) 



61n(l/r2A2 CD ) 12iV c 

and where one would need an appropriate prescription to freeze it when r starts to approach 
1/Aqcd- Now the energy dependence of the saturation momentum reads (Ml EZ1 E3 EH] 



dlnQ 



dY 



xds) i 161 


\x"(ls)] 


1/3 


"x(7s)" 


1/6 1 


2bj s Y 8 








y5/6 



(3-11) 



with £i = —2.33 the leftmost zero of the Airy function. In order to find the behavior of the 



S'-matrix at saturation, we make in Eq. (3.5) the replacement 
1 



a 



Mn(l/^ CD )' (3 ' 12) 

and where clearly this factor has to be moved inside the integrand. Notice that such a replacement 
can be obtained from a "smallest dipole" prescription, which is a natural on^J given the dipole 
splitting (xi,X2) — (xi, z) and (z, X2) we let the coupling run according to a{r 
size of the smallest of the three dipoles involved in the process. We find 



mm ,, with r min the 



d(S) 



Y 



dY 



-I in 
b 



ln(Q^ CD ) 
ln(l/r^ CD ) 



and after the integration by keeping only the most dominant term in Y we obtain 
1 



(S)y ~ exp 



2b 



YlnY 



(3.13) 



(3.14) 



3 In the subsequent sections we shall extend our study to include two more prescriptions: the Balitsky prescription 
[65 which has dominated the phenomenology in the recent years, and a particular daughter dipole prescription 



which is used when one solves the Langevin form of the JIMWLK equation [401 155] . 
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FIG. 5. Numerical extraction of the coefficient a\ in the asymptotic expansion — ln(5(r))y = a{Y hiY + 
a^Y + as In Y + at saturation and for running coupling. 

In the numerical solution shown in Fig. [5] we see that one indeed recovers this dependence in Y 
without the need to evolve too much. However the value of the coefficient differs around 30% from 
its asymptotic value 1/26. The approach to the latter is extremely slow and takes place only at 
tremendously high values of Y. 



IV. THE GAUSSIAN APPROXIMATION 

As said at the end of Sect. [TTJ one method to calculate multi-gluon correlators is to solve the 
Langevin form of the JIMWLK equation on a lattice [39J. The implementation of this approach 
[iO] lead to the remarkable finding that the numerical data for the quadrupole are well-described by 
a Gaussian approximation, more precisely by the extension of the MV model to arbitrary rapidity 
Y. We shall not review here the MV model, but is suffices for our purposes to recall that it is the 
typical initial condition at some Yq. It is a Gaussian functional in the color sources p a and therefore 
also in the field a a in the covariant gauge, but with a kernel which is independent of Y , since it 
refers to a fixed, initial, value of rapidity. Furthermore, it is clear that any Gaussian probability 
distribution involves only a single kernel, thus all high-point correlators can be expressed in terms 
of the 2-point one, e.g. the dipole 5-matrix (S)y, a property which is true in the MV-model. 
Nevertheless, there was no a priori reason for this to happen for the JIMWLK Hamiltonian which 



is highly non-linear due to the Wilson lines in Eq. (2.2) and which arise, as said, from the scattering 



of the emitted gluon with the exisiting background target field. Perhaps a hint was given some time 
ago in [69], where a "random phase approximation" to JIMWLK lead to a Gaussian Hamlitonian, 
even though no explicit reference to high-point correlators was made there. A Gaussian ansatz has 
also been used in [70] in order to estimate 1/N% corrections to the BK equation. 

Still, despite the aforementioned nonlinearities, one can prove that a Gaussian approximation is 
a quasi-exact solution to the JIMWLK equation j3U[32]. The non-linearity is indeed present in the 
dipole equation given in Eq. ( |2.12 ), but not in the quadrupole one, Eq. (2.14), which is linear in 



Q (at the operator level and/or at large- N c ). This is an indication that a Gaussian approximation 
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may be a possible solution with all the non-linearities absorbed in its kernel or, equivalently, in the 
2-point function. Let us shortly review and explain how this happens \%T\ W2\ . 

At saturation, where by definition the target is dense, real emissions are suppressed and the 
virtual part i7 v irt of the Hamiltonian dominates. This part has obviously a Gaussian form, including 
the second term; those adjoint Wilson lines simply transform the "left" functional derivatives to 
"right" ones which act on the lower and upper end-points of the Wilson lines W and V respectively, 
namely 

<%r± vi = (vir^ « = ^ « = w- vie, ^ 

and similarly for the action on V x . 

Now we proceed as in the case of the BK equation. Since the integrand of i? v irt is ^-independent, 
we can integrate the dipole kernel over z in the region 1/Q 2 <C \u — z\ 2 ,\v — z\ 2 <C \u — v\ 2 . 
Again the lower limit is imposed by our approximation, while the upper one is chosen to give the 
dominant logarithmic contribution which is 2 In [(u — v) 2 Q 2 ] . So far our approach would lead to 
a Hamiltonian valid only at saturation, but recalling Eq. ( |3.5| ) we see that this logarithm can be 
expressed in terms of the logarithmic derivative of the dipole w.r.t. Y and with such a replacement 
we finally arrive at 



d]n(Suv /Y 



— I 

] 2 Cf Jui 



A„1ri / jv ( si„.a si~.a z^.a z^.a ) ' (4.2 



This is a Gaussian Hamiltonian which is correct, for finite- Af c , at saturation by construction and 
in the dilute limit as can be inspected. The kernel, which has absorbed the non-linearities, is In- 
dependent (in contrast to the MV model one) and is most easily determined from the BK equation. 
Indeed one can verify that in the Gaussian approximation, and in an arbitrary representation i?, 
we have 

ln(5£,} y = ^ln<£l K ) y , (4.3) 



and insert it for R = F in Eq. (4.2). Therefore, by solving the BK equation, a closed equation 



at large- N c , we can first find the dipole S'-matrix at finite- N c using Eq. (4.3). Then, using Hq 
one constructs evolution equations for high-point correlators which have the advantage to be local 
in the transverse plane, which simply means they are ordinary differential equations in Y with 
Independent coefficients. For example, the quadrupole at large- N C1 on which we will focus in 



Sect.[VlJ reads gUEE] 

(QX1X2X3X4 )Yq 



{Qx 1 x 2 x 3 x 4 )y — V (S Xi x 2 )y (Sx3X 2 )y (S X3 X4)y (S Xi x4)y 



+ i r Y dy 



(Sx 1 x 2 )y (Sx 3 x 2 )y (S X3X4 }y (S XiX4 )y 

(S x lX 3 )y{Sx2X4)y 9 (S XlX2 )y(S X3X4 )y + ( S Xl £C 4 ) y ( ^X 3 X 2 ) V 



. (4.4) 

The quadrupole above obeys the "mirror symmetry" (Q XlX2X3X4 ) Y = (Qx 1 x4x 3 x 2 ) Y i a P ro P ert y 
which in fact holds at finite- N c and beyond the Gaussian approximation and is a consequence of 
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symmetry under time-reversal where time is represented by x~ [42J. This symmetry is preserved by 
the JIMWLK equation due to the two types, left and right, of functional derivatives and suggests 
that the hadron expands symmetrically in the x~ direction during the evolution. At the level of 
the Gaussian approximation only, the quadrupole is also symmetric under the charge conjugation 
(Qx 1 x 2 x 3 x4)y ~ (Qx 2 x 3 x4x 1 )y- ft i s a l so important to check the weak scattering limit of Eq. ( |4.4[ ); 



when all (T XiXj )Y = 1 — (S XiXj )Y are small, then the expansion of Eq. (4.4) to first order in (T)y 
reads 

(Qxi X2X3X4 }y = 1 — (T XiX2 )y + (T XiX3 )y - (T XiX4 )y - (T X2X3 )y + (T X2 x 4 )y - (T X3X4 )y: (4.5) 

which is assuredly the same relation one finds by expanding the Wilson lines of the dipole and the 
quadrupole to order {go) 2 . Let us also point out that a general initial condition where (Q)y is n °t 
determined by (S)y like, for example, the one in [71J, can be accommodated for the quadrupole 



in Eq. (4.4). 



In some scenarios one can perform analytically the ^-integration in an equation like Eq. (4.4) 
and therefore arrive at a functional expression for a high-point correlator in terms of the dipole 
which is local in Y . Such a situation is realized for all the simple configurations studied in 



and independently of whether the coupling is running or not. Furthermore, using a "separability" 



property of the Gaussian kernel jH] S2] in Eq. (4.2) one can show that this is also the case in 



fixed coupling evolution for an arbitrary configuration. For example, the quadrupole in Eq. (4.4) 
assuming MV model initial conditions at Yq reads 

/A V _ m [(Sx 1 X2)y(Sx3X4)y/(SxiX3)y(Sx2X4:)y\ /q \ /C \ 

\Wx 1 X2X 3 X4/y - , r/A v /A \ //A \ /A \ ] \^X!X 2 / Y \^3^ /Y 

111 [\^XlX2 /Y \^X 3 X4/Y I \^XiX4/y\ Xc 2X3/y\ 

In [(«S f agias 4 )y(^ag2a5 3 )y/ (^iz 3 )y {$x 2 X4 )y] / A \ /A \ (A a \ 

i r/A \ /A \ //A \ /A \ ^ \^x lX4 /Y\^X2X3/y V 4 - D J 

m l\dxiX4 /Y\ X 2X3 /y/ \^XlX2 /y\^X3X4 /yj 



+ 



At the formal level this expression was first derived in the MV model [T7] at large- N c and later 
on generalized at finite- N c [33J. One non-trivial achievement in [41] [42] was to show when and 
why it remains (approximately) valid after quantum evolution has been taken into account. Like 



Eq. (4.4) it also reduces to Eq. (|4.5|) in the weak scattering limit. We note that for an arbitrary 



configuration in running coupling evolution one is supposed to use Eq. (4.4), since the arguments 
leading to Eq. ( |4.6[ ) do not go through (except in the case that the scale of the running coupling 
is taken to be Q s ). Still we shall see that, so long as rQ s is not becoming too large, one can also 



use the simpler version Eq. (4.6) for practical purposes. 



At finite- N c there is operator mixing and the analogue of Eq. (4.6) involves the diagonalization 
of a matrix. For example, the quadrupole mixes with the two dipoles operator, while the phe- 
nomenologically interesting 6-point operator Qx 1 x 2 x 3 x4Sx4x 1 5 appearing in the quark-gluon double 
inclusive cross section, mixes with two more operators and one can show that the emerging cubic 
equation has a unique real solution. It goes without saying, analytic expressions for high-point cor- 
relators are invaluable in order to reduce the numerical cost for Fourier transforming to momentum 
space and obtain the desired cross sections. 
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FIG. 6. Precise solution to the BK and MFA equations for fixed coupling. 

V. A MEAN FIELD EQUATION 

One way to test the accuracy of the approximation scheme is by requiring that the Gaussian 



Hamiltonian Hq in Eq. (4.2) coincides with the "averaged" JIMWLK Hamiltonian H in Eq. (2.2) 
[231 H2]. Then it is an easy and straightforward exercise to obtain a closed non-linear equation for 
the dipole in the adjoint representation A which reads [12] 



9 ( S X 1X2 )y _ a f KA ( S X 1X2 ) 



- [ M \ b viK2/Y (cA _j_ a A _ cA _-,\ 

QY n JVl XlX 2 Z . , qA \ ^ XlZ ZX2 X 1 X 2 L /Y' 

J Z \ ^XlX 2 /Y 



(5.1) 



We shall study the solution to the above in the large- 7V C limit, in which (S£ iX2 )y = (S Xi x 2 )y an( i 
then compare to that of the BK equation. A quick inspection shows that Eq. (|5.1[) reduces to 



the BFKL equation in the limit of weak scattering and to Eq. (3.5), like the BK equation, in the 
regime of strong scattering. 

This expected analytic behavior is in fact observed when we compare the numerical simulations 
of the BK equation Eq. (2.12) to the simulations of the MFA equation (5.1). In Fig. [6] we do see 



that the two equations, using the same initial condition, agree in the regime of weak scattering. 
When approaching the saturation regime the two solutions start to separate and, even though 
(S)y K and (5)y FA are not the same at saturation, their logarithms agree quite well as we can 
easily infer from the right plot in the upper panel and the left plot in the lower panel of Fig. [6] and 
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FIG. 7. Precise solution to the BK and MFA equations for running coupling. Above: smallest dipole 
prescription. Below: Linear scale plots for various schemes for the running of the coupling. 



as was advocated at the beginning. 

Eq. (5.1) is equally good in running coupling evolution so long as we adopt, as we should, the 
same prescription for the running of the coupling in both the BK and JIMWLK equations and 
indeed this is what is shown in Fig. [7[ In the upper panel we show the results for the smallest 
dipole prescription, introduced earlier in Sect. Ill, which are similar to the fixed coupling ones 
in the upper panel of Fig. Eq. Q. Even though not shown in a logarithmic scale, to avoid a 
proliferation of plots, the results for the daughter dipole prescription as in [101 [55] and the Balitsky 
prescription as in [65], are almost identical in shap^j In fact this property can be inferred from 
the linear plots in the lower panel of Fig. [7[ 

We shall see in the next Section that, when we express higher-point correlators in terms of the 
dipole, the accuracy is not restricted to the logarithmic level. 



VI. QUADRUPOLE CONFIGURATIONS 



Ideally, one would like to solve the quadrupole equation in general. In principle this is possible 
since it is a closed equation (at large- iV c ), but one understands that in practice it is rather com- 

4 Here we mention the observation of a "peculiar", and perhaps unphysical, feature already pointed out in [72]: the 
solutions (for both BK and MFA) with the Balitsky prescription evolve slower than those with the smallest dipole 
one. In 65, 66 it was shown that the Balitsky prescription is equivalent to the smallest dipole one in the limits 
where the dipoles sizes are very different. However one can see that when the daughter dipoles are large, that is 
when \xi — z\ c± \x2 — z\ ^> \x\ — x^\ then one finds aeai — ol{\x\ — a?2 1)(1 — 4a(|sci — z\)) and therefore this 
convergence is very slow. 
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(a) (b) 

FIG. 8. (a) The "line" configuration of the quadrupole for which 7*12 = 7*14 = 7*23 — 7*34 with 7*13 = r24 = 0. 
(b) A triangle configuration emerging from the evolution of the "line" . 



plicated, because of the large number of transverse variables (8 in general, 6 if we assume impact- 
parameter independence and 5 if we also impose rotational invariance) on which the quadrupole 
depends and because of the non-locality in the transverse plane. Therefore we have to be more 
modest in our goals, and what we shall do is to average the quadrupole with a Gaussian wave- 
function, thus obtaining an evolution equation for the dipole. The solution to the latter will be 
compared with the solution to the BK equation, and such a comparison should give us a good 
estimate of the validity of the Gaussian approximation. Still, many coordinates appear in such a 
general equation and therefore we shall consider only a couple of special configurations. We shall 
immediately see that the one of them can be analytically investigated. 

a. The "line" configuration To start, let us consider the "line" configuration of the quadrupole, 
as shown in Fig. |8j(a), where we put the two quarks at the same coordinate and similalry for the 
two antiquarks. This is the simplest possible one, since it is characterized by only one non-vanishing 



distance, and its evolution, according to Eq. (2.14), is given by 



d(Q XlX2XlX2 )Y _ a f M /A A A A A _ A2 \ / fi x 

QY — ^ J JVl XlX2Z\°XlZ Kc 6 ZX2X1X2 ^ °ZX2 Kc CXlZXlX2 K °CX\X2X\X2 ° X\X2 I Y ' \ ) 

Let us assume the large- N c limit, so that we can factorize (SQ)y —> (S)y(Q)y an d (S 2 )y —> (S)y- 
Then the above equation becomes a closed equation for the quadrupole, where the dipole is known 
from the solution to the BK equation. Still, as expected, this is not a closed equation for the 
particular quadrupole configuration, since a more general one appears in the real terms on the 
r.h.s. of the equation. It has the shape of a triangle as shown in Fig. |8j(b); only the two quarks are 
at the same point with the two antiquarks being separated, or vice versa. As we have said above, 



we shall assume a Gaussian average over Eq. (6.1), and use (4.6) in order to obtain an equation 



for the dipole. To this end we need the corresponding average for the configurations appearing in 



Eq. (6.1) and which are given by or obtained from (always in the large- N c limit) 



{Qx 1 zx 1 x 2 )y = {S XiX2 )y{S XiZ )y [l + ^((S xiX2 )y(S XiZ )y/ {S X2Z )y)] , (6.2) 



and then Eq. ( |6.1[ ) leads to 
d(S XlX2 ) Y _ ol 



dY 



which is nothing else than the BK equation. We would also like to stress here, without going 
through the details of the derivation, that the above is still true even when N c remains finite; it is 
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(a) 



(b) 



FIG. 9. (a) The "square" configuration of the quadrupole for which 7*13 = r u and r 2 3 = r 2 4- The average 
value of Q depends only on the distances depicted by straight lines, (b) A configuration emerging from the 
evolution of the "square" . 



a straightforward exercise to show that taking the Gaussian average of Eq. ( |6.1[ ) we are lead to 
d\n(S XlX2 ) Y 



dY 



1 a r 



M 



{S Xi z}y{Szx 2 )y 

(S Xi x 2 )y 



(6.4) 



which is the Gaussian average of the dipole equation at finite- N c . Moreover, let us note that one 
arrives again at the above when considering the Gaussian average in the evolution equation of 

(s < x 1 x 2 )y- 

All this is a strong indication that the Gaussian approximation should be a very good approx- 
imation at least to the particular configurations under consideration. However, one should not 
draw the conclusion that it is the exact answer. Indeed, one can proceed to assume a Gaussian in 
the evolution equation of the triangle quadrupole configuration in Fig. |8j(b) to find that this time 
it does not reduce to the BK equation for the dipole. 

b. A "square" configuration Now let us study the configuration shown in[9j(a) where the four 
fermions are located at the corners of a squar^J Nevertheless, we prefer to draw the diagonals, 
for reasons that we immediately explain. Clearly, a "square" configuration is simple enough in the 
sense that there are only two different distances (the side and the diagonal) between the fermions. 
This is one of the special class of configurations where the quadrupole can be written as a product 
of two dipoles in the Gaussian approximation JJTJH2]. More precisely, for any value of N c , one has 



(Qx 1 x2x 3 x 4 )y = (S XiX2 )y(S X3XA )y = (S(R))y, (6.5) 

where with R we denote the length of the diagonal. The r.h.s. of the quadrupole equation involves 
quadrupole configurations more complicated than the "square", e.g. {Q Z x 2 x3xa)y ^ where z is lo- 
cated anywhere in the 2-dim transverse plane as shown in Fig. [9j(b). This also means that our 
test of the Gaussian approximation is probing not only the simple square configuration, but also 
the wider sample of configurations represented by the one in Fig. [9j(b). Using the symmetries of 
the configuration under consideration we can regroup the terms on the r.h.s. of the equation which 



5 Notice that this is different than the "square" configuration studied in [40]. There, the two quarks are at the edges 
of the one diagonal and similarly for the antiquarks. Here, they are at the edges of the same side. 
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FIG. 10. Solution to the BK and the reduced quadrupole equation for the configuration in Fig. [9] for fixed 
coupling. 



a 



1 



simplifies to 

d(S(R)) Y 

OY 2tt (S(R)) y 

In general, the quadrupole in the Gaussian approximation and at large-iV c is already given in 



L {Mxi 



X2Z +M X -M 



] [{S XiZ )y (Q ZX2X3X4 )y-(S(R)) 2 y }. (6.6) 



Eq. (4.4) or its simplified version Eq. (|4.6|) and we shall use the latter to rewrite (6.6) as a closed 



equation for (S(R))y- Then considering {Qzx 2 x^xa)y in Eq. (6.6), some minor cancelations occur, 
since (S X2X3 )y = (S X2X4 )y, but one cannot go far in simplifying the general expression in Eq. (4.6). 



As a check of our manipulations let us note that in the weak scattering limit Eq. (6.6) reduces to 
the BFKL equation for (S(R))y — (S XiX2 )y, while in the strong scattering regime one recovers to 



Eq. (3.5). All this is natural if the Gaussian Hamiltonian is expected to be a good approximation 
scheme. 



The numerical results for fixed coupling evolution are shown in Fig. 10 and they should be 
compared to those of the MFA in Sect. |V| shown in Fig. [6) It is clear now that the curves arising 



from the solution to the BK equation and Eq. (6.6) for the square configuration almost fall on top 



of each other in all kinematic regimes. This is a strong indication that the Gaussian approximation, 
and more precisely the extrapolation to arbitrary Y of the MV model, is a quasi-exact solution to 
the JIMWLK equation for fixed coupling evolution. 

The situation is almost the same for running coupling evolution and for all the prescriptions 



used, which are those adopted in Sect. N] This is exhibited in Fig. 11 where one sees that in 
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FIG. 11. Solution to the BK and the reduced quadrupole equation for our particular configuration for 
running coupling. 



practice Eq. (4.6) provides an excellent approximation to the quadrupole. Some deviations start 
to occur when rQ s becomes large and in fact we have seen in the numerical solutions that these 
deviations grow and extend deeper in the saturation region when rapidity the Y increases, although 
an agreement similar to the one observed at fixed coupling is eventually recovered deep in the 



saturation region. Most likely this discrepancy is to be attributed to the use of Eq. (4.6) instead 
of the proper expression for running coupling evolution given in Eq. ( |4.4| ) (cf. the discussion above 
and below Eq. ( |4.6| )). This is also supported by the fact that the deviations are stronger (milder) 
for the Balitsky (daughter dipole) prescription as a consequence of the stronger (milder) variation 
of the coupling with the scale. Thus, it should be interesting to check if the discrepancy goes away 



if one makes use of Eq. (4.4), even though this looks to be mostly an academic problem at least 



for this particular configuration. 



VII. CORRELATORS WITH OPEN COLOR INDICES 

So far we have dealt with observables which involve only traces over products of Wilson lines. 
However, there are quantities which require the knowledge of correlators with open color indices, 
like the energy density and its fluctuations immediately after the collision of two heavy ions, and 
we shall present here a way to calculate them. We start from the simplest possible case which 
is the operator {Vx 1 V X2 ) i ^ with i,j = 1, 2, 7V C . Its average value is proportional to the color 
structure 5^, so that we must have 

(«y X2 ) iJ ) Y = 5 iA^2) Y - (7.i) 

Indeed, if we contract the above with 5i« we arrive at the average of the definition of the dipole 



S Xl x 2 as given in Eq. (2.7). This may look rather trivial for the moment but we shall see that the 



method works for less simple operators with more open color indices. It is also important to point 



out that the color structure in Eq. (7.1) is preserved under evolution, either JIMWLK or Gaussian. 



At this level, one needs to be careful and use the more general form of the JIMWLK Hamiltonian 

H = s* L*- 4 ( x + ^ - m - w-T w <72) 
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with the kernel 



(u 



(V 



(u — z) 2 (v — z) 2 ' 



(7.3) 



when acting on Wilson lines with open indices. Considering, as an example, only the last (the one 
involving VzV v ) of the four terms in the parenthesis of Eq. (7.2), it is a straightforward exercise 
to show that 



[H{VlV a 



x 2 Jij/Y 



4th 47T ' 



VX 2j 



5ux! (^x 1 z{Vz^x 2 ) ij)y ^ux 2 ((^x 1 ^z) ij^zx 2 ) 



Y 



ij/Y 



(7.4) 



Now we employ Eq. (7.1) in the last term in the bracket, while we furthermore assume for the first 
term that 



(^x 1 z(Vz^x 2 ) ij)y ~ 3ij(Sx 1 zSzx 2 )Yi 

and similarly for the second one. Making use of 



/ 



fcuvz {$ux\ 



J UX 2 ) \ u VXl 



J vx 2 ) 



M X1 



X 2 Zi 



(7.5) 



(7.6) 



l uv 

we can write 



M X1 



= Sij(HS XlX2 ) Y . (7.7) 

1 4th 



We can follow the same procedure for all terms of the JIMWLK Hamiltonian to finally arrive at 
(H(VIV X2 ) ) Y = (HSij S XlX2 ) y , (7.8 



which was our original claim. It means that if Eq. (7.1) and equations of the form (7.5) are valid 
at rapidity Y, then Eq. (7.1) will hold at Y + Ay. Subsequent action of the Hamiltonian will 
generate operators with more and more Wilson lines, but the only structures that can emerge are 
SijO (such a term in the above example is not shown in Eq. (7.4) but is generated by the first term 
of the Hamiltonian) or (VW . . . VW)ijO, where 6 is an operator with no open color indices. This 
will happen because we are interested in Wilson lines in the fundamental representation; the two 
functional derivatives of the Hamiltonian give rise to the color structure {t a )ki{t a )mn which can be 
expressed in terms of Kronecker deltas, through the Fierz identity. Thus one finally arrives at the 
conclusion that (with S^ 2n ^ given in Eq. (|2.9|)) 



Y 



(7.9) 



will be true for any O with no open color indices and at any value of the rapidity Y, since it is 
trivially valid for vanishing color gauge field. 

Let us now go to a more interesting case by considering the operator (V Xl V X2 ) .^(V^^ 3 V X4 ) kl O. 
Similar reasoning leads us to expect the color structure 



((K^L {VlV X4 ) kl 6) Y = (A) Y + 6 u 6 jk {b) y + s lk s Jt {c) 



(7.10) 
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and by contracting with 5ij8kh $U&jk an d SikSji we arrive at the 3x3 inhomogeneous system of 
linear equations 



(7.11) 



where we have defined the operator R XlX2X3X4 = (l/N c ) (V^ V^)^ (V^ V X4 ) ir Solving Eq. ( |7. 



X i r 










1 N c 1 




{B)y 






1 1 JV C _ 




_{C)y_ 




(-RX1X2X3X4&") y 



we obtain the gauge-invariant operators 

N C {N C -\~ '\-)( k SxiX2^X3X4^')y ~ (QxiX2X3X4@ S )y ~ 



{B)y 

(C)y 



X2X3xJ3/Y 



Y (N C + 2)(N C -1) 

(N c + 1) {^Qxi X2X3X4&) y ~ ^\^X\X2X3X4&) y ~ N c\Sx\X2^X3X4^) y 

~ (iV c + 2)(7V c -l) : 

(N C + 1)(R 

x\X2X3xJ^)y i^Q X1X2X3X 4^^) y N c (S XiX2 S X3X4 (D/y 



(7.12) 
(7.13) 
(7.14) 



(N C + 2)(N C -1) 

This is general and one would need to calculate separately the expectation values of the operators 
SSO, QO and RO. Of course considerable simplifications are expected to occur in the Gaussian 
approximation; indeed, one finds that SSO, QO and RO are such that (C) y in Eq. (7.14) vanishes 
and, thus, the color structure S^Sji does not appear any more in Eq. (7.10). Furthermore, by 
eliminating RO, we can simplify (A) y and (B} y to finally arrive at 

// T/ t T/ \ / T /t t/ \ A\ r r N c(^x 1 x 2 S X 3X4^)y ~ (Q^iX2X3X4^) Y 
\\ V k V X2)ij\ V x\y X 4)kl /Y = d H d kl N 2 _ x 

N c{Q x \X2X3X 4^) y ~ N c\Sx\X2 &X3X4 



+ 8u8jk 



N*-l 



(7.15) 



In the large- N c limit one can neglect the term involving the quadrupole in the first fraction and 
also set N% — 1 « N% in the denominators. However, one should be careful not to perform any 
further large- N c approximation at this stage; a measurable quantity will involve contractions over 
the color indices in Eq. (7.15) and such a procedure can alter the N c counting. 



VIII. CONCLUSIONS AND PERSPECTIVES 



Analytic expressions for multi-gluon correlators in the high energy limit are indispensable in 
order to reduce the numerical cost for obtaining the cross sections related to multi-particle pro- 
duction. Such expressions arise in the context of the Gaussian approximation to the JIMWLK 
evolution equation, and in this work we have studied the validity of such an approximation. Our 
results are already succinctly expressed in the last paragraph of the Introduction (Sect.]!]). In short, 
we have confirmed that the Gaussian approximation provides a quasi-exact solution and in partic- 
ular we have shown that this is the case independent of the prescription used to set the scale in 
the argument of the running coupling. 

Still, there remain a few things which could be addressed in a future work. It is important 
to explore how solid the Gaussian approximation is when we use more general initial conditions, 
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like, for example, those in [71J and we do believe that the answer to this question is positive. 
Furthermore, as discussed in Sect.|VI| it would be interesting to check if the use of the more general 



expressions in the Gaussian approximation, like the one in Eq. (4.4), instead of the extrapolation of 



the MV model to arbitrary Y, like in Eq. (4.6), improve the accuracy in running coupling evolution. 
Last, but not least, one should perhaps examine how and if the approximation applies to the case 
of unequal rapidity correlations. 
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Appendix: Numerical implementation 

In this Appendix, we sketch the main points of the numerical techniques used to solve the evolu- 
tion equations in this paper. Most of them are common to the BK equation, that is the factorized 



version of Eq. (2.12), the MFA equation Eq. (5.1) and the equation of the square configuration 



Eq. (6.6), so we will focus on the BK equation below and comment on the other cases after. 

For the integration of the right-hand side of the equation, the first step is to use the mirror 
symmetry along the bisector of the external dipole line xi, x 2 and limit the integration to the 
half-plane containing x\. In that region, we integrate in polar coordinates around x\ and easily 
get 

/ m x1X2Z = r di r~ 9o{£) ( * 1 = x f 2 (a.i) 

where £ — ln(l/(xi — z) 2 ) and the bounds on the 9 integration — with 9 measured from the 12 
axis — depend on £. The radial and angular integrations are then discretised. We typically use 
16 points in 9 and 1440 points in £ with —40 < £ < 140. We have checked that for all the plots 
presented in this paper, the discretisation in £ was fine enough to reach the quoted accuracy. 

Then, since deep at saturation (5(r))y becomes extremely small when rapidity increases, we 
decided to work instead with s(r, Y) = ln(l/(*S(r))y) and all the equations are rewritten in terms of 
s. To achieve maximal numerical precision, terms that cancel the logarithmic ultraviolet divergence 
when z approaches x\ are grouped. 

Knowing the right-hand side of the equation, the numerical evolution in rapidity is then handled 
using a standard fourth-order Runge-Kutta method. 

This method applies straightforwardly to the BK and MFA equations. For the equation of the 
square configuration, the bisector of the (xi, x 2 ) dipole goes through X4 and one may wonder about 
the possible ultraviolet divergence when z approaches X4. It is easy to notice that this appears 
under the form of the combination of BFKL kernels A4 XlX4Z — M X2X4Z i which is only linearly 
divergent for z #4, multiplied by a combination of the dipole amplitude S(z) that vanishes 
when z X4, giving an overall smooth behaviour. 
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Finally, let us notice that the formulation of the evolution in terms of s(r,Y) = ln(l/(S r (r))y) 
also has its limitations if Y becomes too large. This is because, in order to obtain the full asymp- 
totic geometric scaling behaviour, one needs to describe correctly the dilute tail over the whole 
geometric scaling region. Since this extends to smaller and smaller (T(r))y, with (T(r))y « 
s(r, y), when rapidity increases, a better precision in the dilute region is obtained by working 
with ln(l/(T(r))y). This can be combined with the dense region by using a mixed variable 
X(r, Y) = max{ln(2/(S r (r))y), ln(2/(T(r))y)}. This can be implemented using the techniques 
above and keeping track of the threshold £q where the transition between X = ln(2/(*S(r))y) and 
X = ln(2/(f (r))y) occurs, i.e. where (S(£ ))y = (T(£o))y = 1/2. It is this method that we have 
developed to verify the Levin- Tuchin formula. 



[1] BRAHMS Collaboration, I. Arsene et al, Phys.Rev.Lett. 93, 242303 (2004), nucl-ex/040 3005 . 

[2] STAR Collaboration, J. Adams et al, Phys.Rev.Lett. 97, 152302 (2006), |nucl-ex/0602011| . 

[3] STAR Collaboration, E. Braidot, Nucl.Phys. A854, 168 (2011), [1008.3989]. 

[4] PHENIX Collaboration, A. Adare et al, Phys.Rev.Lett. 107, 172301 (2011), [1105.5112]. 

[5] D. Kharzeev, Y. V. Kovchegov and K. Tuchin, Phys.Rev. D68, 094013 (2003), |hep- ph/0307037l . 

[6] J. L. Albacete, N. Armesto, A. Kovner, C. A. Salgado and U. A. Wiedemann, Phys.Rev.Lett. 92, 
082001 (2004), [hep-ph/0307179] . 

[7] J. P. Blaizot, F. Gelis and R. Venugopalan, Nucl.Phys. A743, 13 (2004), |hep-ph/0402256l . 

[8] E. Iancu, K. Itakura and D. Triantafyllopoulos, Nucl.Phys. A742, 182 (2004), [ hep-ph/0403103| . 

[9] J. L. Albacete and C. Marquet, Phys.Lett. B687, 174 (2010), [1001.1378]. 

[10] T. Altinoluk and A. Kovner, Phys.Rev. D83, 105004 (2011), [1102.5327]. 

[11] J. Jalilian-Marian and A. H. Rezaeian, Phys.Rev. D85, 014017 (2012), [1110.2810]. 

[12] G. A. Chirilli, B.-W. Xiao and F. Yuan, Phys.Rev.Lett. 108, 122301 (2012), [1112.1061]. 

[13] P. Tribedy and R. Venugopalan, Phys.Lett. B710, 125 (2012), [1112.2445]. 

[14] A. Mueller and S. Munier, 1206.1333. 

[15] J. L. Albacete, A. Dumitru, H. Fujii and Y. Nara, Nucl.Phys. A897, 1 (2013), [1209.2001]. 

[16] A. H. Rezaeian, 1210.2385. 

[17] J. Jalilian-Marian and Y. V. Kovchegov, Phys. Rev. D70, 114017 (2004), hep-ph/0405266 . 

[18] C. Marquet, Nucl. Phys. A796, 41 (2007), [0708.0231]. 

[19] K. Tuchin, Nucl.Phys. A846, 83 (2010), [0912.5479]. 

[20] J. L. Albacete and C. Marquet, Phys. Rev. Lett. 105, 162301 (2010), [1005.4065]. 

[21] A. Stasto, B.-W. Xiao and F. Yuan, Phys.Lett. B716, 430 (2012), [1109.1817]. 

[22] T. Lappi and H. Mantysaari, 1209.2853. 

[23] E. Iancu, A. Leonidov and L. McLerran, |hep-ph/0202270| 

[24] Y. V. Kovchegov and K. Tuchin, Phys.Rev. D65, 074026 (2002), [hep-ph/01 11362] , 

[25] A. Kovner and M. Lublinsky, JHEP 0611, 083 (2006), hep-ph/0609227 . 

[26] F. Gelis, T. Lappi and R. Venugopalan, Phys.Rev. D79, 094017 (2009), [0810.4829]. 

[27] K. Dusling, F. Gelis, T. Lappi and R. Venugopalan, Nucl.Phys. A836, 159 (2010), [0911.2720]. 

[28] Y. V. Kovchegov and D. E. Wertepny, 1212.1195. 

[29] F. Dominguez, C. Marquet, A. M. Stasto and B.-W. Xiao, 1210.1141. 

[30] L. D. McLerran and R. Venugopalan, Phys. Rev. D49, 2233 (1994), hep-ph/9309289 . 

[31] L. D. McLerran and R. Venugopalan, Phys. Rev. D49, 3352 (1994), |hep-ph/9311205l . 



25 



J. P. Blaizot, F. Gelis and R. Venugopalan, Nucl. Phys. A743, 57 (2004), hep-ph/0402257 . 

F. Dominguez, C. Marquet, B.-W. Xiao and F. Yuan, Phys. Rev. D83, 105005 (2011), [1101.0715]. 

J. Jalilian-Marian, A. Kovner, A. Leonidov and H. Weigert, Nucl. Phys. B504, 415 (1997), hep- 

lph/9701284) . 

J. Jalilian-Marian, A. Kovner, A. Leonidov and H. Weigert, Phys. Rev. D59, 014014 (1998), hep- 
|ph/9706377l - 

H. Weigert, Nucl. Phys. A703, 823 (2002), [hep-ph/0004044l . 



E. Iancu, A. Leonidov and L. D. McLerran, Nucl. Phys. A692, 583 (2001), hep-ph/0011241|. 

E. Iancu, A. Leonidov and L. D. McLerran, Phys. Lett. B510, 133 (2001), [ hep-ph/0102009| . 

J.-P. Blaizot, E. Iancu and H. Weigert, Nucl. Phys. A713, 441 (2003), |hep-ph/0206279| . 

A. Dumitru, J. Jalilian-Marian, T. Lappi, B. Schenke and R. Venugopalan, Phys. Lett. B706, 219 

(2011), [1108.4764]. 

E. Iancu and D. Triantafyllopoulos, JHEP 1111, 105 (2011), [1109.0302]. 

E. Iancu and D. Triantafyllopoulos, JHEP 1204, 025 (2012), [1112.1104]. 

E. Levin and K. Tuchin, Nucl. Phys. B573, 833 (2000), |hep-ph/9908"317] . 

A. H. Mueller and G. Salam, Nucl.Phys. B475, 293 (1996), |hep-ph/9605302] , 

E. Iancu and L. D. McLerran, Phys. Lett. B510, 145 (2001), |hep-ph/0 103032] , 

A. H. Mueller, Nucl. Phys. B643, 501 (2002), h ep-ph/0206~216] . 

E. Iancu and A. H. Mueller, Nucl. Phys. A730, 494 (2004), |hep-ph/0309276l . 

I. Balitsky, Nucl. Phys. B463, 99 (1996), |hep-ph/9509348] . 

Y. V. Kovchegov, Phys. Rev. D60, 034008 (1999), [hep-ph/990128l] . 

Y. Hatta, E. Iancu, K. Itakura and L. McLerran, Nucl. Phys. A760, 172 (2005), |hep-ph/0501171| . 
A. H. Mueller, Nucl.Phys. B415, 373 (1994). 



I. Balitsky, hep-ph/0101042. 



D. N. Triantafyllopoulos, Acta Phys. Polon. B36, 3593 (2005), hep-ph/0511226| . 
K. Rummukainen and H. Weigert, Nucl. Phys. A739, 183 (2004), hep-ph/0309306| . 
T. Lappi, Phys.Lett. B703, 325 (2011), [1105.5511]. 

A. H. Mueller and D. N. Triantafyllopoulos, Nucl. Phys. B640, 331 (2002), hep-ph/0205167| . 
S. Munier and R. B. Peschanski, Phys.Rev. D69, 034008 (2004), |hep-ph/0310357l . 
L. V. Gribov, E. M. Levin and M. G. Ryskin, Phys. Rept. 100, 1 (1983). 

E. Kuraev, L. Lipatov and V. S. Fadin, Sov.Phys.JETP 45, 199 (1977). 
I. Balitsky and L. Lipatov, Sov.J.Nucl.Phys. 28, 822 (1978). 

A. M. Stasto, K. J. Golec-Biernat and J. Kwiecinski, Phys. Rev. Lett. 86, 596 (2001), hep-ph/0007192| . 

E. Iancu, K. Itakura and L. McLerran, Nucl. Phys. A708, 327 (2002), [hep-ph/0203137| . 

S. Munier and R. B. Peschanski, Phys. Rev. Lett. 91, 232001 (2003), [hep-ph/0309177| . 

Y. V. Kovchegov and H. Weigert, Nucl.Phys. A784, 188 (2007), |hep-ph/0609090] . 

I. Balitsky, Phys.Rev. D75, 014001 (2007), |hep-ph/ 06 09105] , 

I. Balitsky and G. A. Chirilli, Phys.Rev. D77, 014019 (2008), [0710.4330]. 

D. N. Triantafyllopoulos, Nucl. Phys. B648, 293 (2003), |hep-ph/0209~12l] . 
G. Beuf, 1008.0498. 

E. Iancu, K. Itakura and L. McLerran, Nucl. Phys. A724, 181 (2003), hep-ph/0212123| . 

Y. V. Kovchegov, J. Kuokkanen, K. Rummukainen and H. Weigert, Nucl. Phys. A823, 47 (2009), 
[0812.3238]. 

A. Dumitru and E. Petreska, Nucl.Phys. A879, 59 (2012), [1112.4760]. 

J. Berger and A. M. Stasto, Phys.Rev. D84, 094022 (2011), [1106.5740]. 

D. Binosi and L. Theussl, Comput. Phys. Commun. 161, 76 (2004), |[hep-ph/0309015| . 

D. Binosi, J. Collins, C. Kaufhold and L. Theussl, Comput. Phys. Commun. 180, 1709 (2009), 



26 



[0811.4113]. 



